Umbilical points of a non-Gaussian random field 



A.M. Turner, 1 T.H. Beuman, 2 and V. Vitelli 2 'E 

'institute for Theoretical Physics, Universiteit van Amsterdam, NL 1090 GL Amsterdam, The Netherlands 
2 Instituut-Lorentz for Theoretical Physics, Leiden University, NL 2333 CA Leiden, The Netherlands 

(Dated: November 6, 2012) 

Random fields in nature often have, to a good approximation, Gaussian characteristics. For such 
fields, the relative densities of umbilical points - topological defects which can be classified into three 
types - have certain fixed values. Phenomena described by nonlinear laws can however give rise to a 
non-Gaussian contribution, causing a deviation from these universal values. We consider a Gaussian 
field with a perturbation added to it, given by a nonlinear function of that field, and calculate 
the change in the relative density of monstars. This allows us not only to detect a perturbation, 
but to determine its size as well. This geometric approach offers an independent way of detecting 
non-Gaussianity, which even works in cases where the field itself cannot be probed directly. 



Random surfaces can be characterized by a variety of 
statistical measures that are geometrical and topologi- 
cal in character. The most well-known are perhaps the 
Minkowski functionals and the statistics of critical points. 
In this paper, we focus on a class of singular points of 
the surface, known as umbilics, that do not depend on 
how the surface is oriented in space. In order to under- 
stand the geometrical meaning of umbilical points, imag- 
ine drawing at every point on the surface the two prin- 
cipal directions, along which its curvature is maximal or 
minimal. At some locations the principal directions can- 
not be defined, because the curvature is the same along 
all directions - these special points are called umbilics. 
As we shall see, umbilical points are topological defects 
with an index of ±1/2. 

This geometrical construct is very useful in a num- 
ber of physical contexts. In statistical optics, the surface 
may represent a curved wavefront that emerges when a 
plane wave is passed through an inhomogeneous refract- 
ing medium. In this mapping, the normals to the surface 
are light rays and the umbilical points correspond to the 
regions where the wave attains its maximal intensity. In 
two-dimensional elasticity or fluid flow, the surface can 
represent a potential function of two variables, whose sec- 
ond derivatives define a shear field that corresponds to 
the principal curvature directions of the surface. The 
points where the shear field vanishes are the umbilical 
points. 

The umbilical points of a surface can be classified into 
three types: lemons, monstars and stars (see figure 1). A 
striking statistical feature of surfaces whose height fluc- 
tuates spatially like an isotropic Gaussian random field is 
that the densities of the three types of umbilics have fixed 
ratios, which are universal numbers [H, 0|- This property 
can therefore be used to test whether a given isotropic 
field is Gaussian; if for a given field h the relative densi- 
ties are found to differ from the universal values, one may 
immediately conclude that the field under consideration 
is not an isotropic Gaussian one. Crucially such a test 
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requires only that the line field corresponding to the prin- 
cipal curvature directions is measurable - the statistics 
of the scalar height field from which the curvature direc- 
tions are derived can be probed without being directly 
observed. 

To give an example of a case where the near-Gaussian 
field of interest is not directly observable, consider the 
phenomenon of weak gravitational lensing Q. As stip- 
ulated by the theory of general relativity, matter bends 
spacetime, which also affects light rays. The light from 
a distant galaxy for instance, does not come to us in a 
straight line, due to the presence of matter between that 
galaxy and us. As a result, we see a distorted image of the 
galaxy. In general, a circular object will look like an el- 
lipse. While most of the matter in the universe is believed 
to be made up of dark matter which we cannot (yet) de- 
tect, the shear field can be detected. The near-Gaussian 
field in this case is obtained by projecting the mass onto 
the sky, along the lines of sight. This is called the pro- 
jected gravitational potential. On large scales, this field 
is approximately Gaussian by virtue of the central limit 
theorem, since the projection involves summing over a 
lot of regions which are randomly distributed. On smaller 
scales however, interactions can give rise to non-Gaussian 
contributions. If we interpret the projected gravitational 
potential as a (near- Gaussian) surface, then the shear 
direction corresponds to the principal direction of this 
surface 4]. At the umbilical points the shear direction 
cannot be defined, hence the amplitude of the shear field 
must vanish - at these locations in the sky, a circular 
light source still appears circular. 

Another example of a physical process in which umbil- 
ical points can prove their usefulness is in the context of 
optical speckle fields. These fields arise for example when 
a coherent beam of light scatters from a rough surface. 
Since the many reflected waves become superimposed, 
this produces a random pattern of intensity with approx- 
imately Gaussian statistics. In this case, it is the points 
of circular polarization that can be identified as umbil- 
ical points. The relative densities of the various types 
of umbilical points have been found to match the theo- 
retical predictions in experiments [5J. A speckle field is 
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not always Gaussian. First, when the surface is not that 
rough, the superposition of the reflected waves will not 
be sufficiently random. Second, a light beam could be 
transmitted through a random medium to map out the 
statistics of its index of refraction 

As a final example, consider a two-dimensional ferro- 
magnet above the critical temperature, for which we mea- 
sure the ^-component of the magnetization as a func- 
tion of x and y. As long as the correlation length of 
the magnetization is smaller than the resolution of the 
measurement, many independent regions are averaged to- 
gether and a Gaussian signal results. However, as one 
approaches the critical point, the correlation length di- 
verges, resulting in non-Gaussianity. 

Other contexts in which umbilical points can offer a 
window for non-Gaussianity include polarization singu- 
larities in the cosmic microwave background [6j-[9( , topo- 
logical defects in a nematic [HI EJ and a superfluid near 

criticality (HE!. 

Testing whether the three types of umbilical points oc- 
cur in their prescribed ratios can thus reveal whether a 
non-Gaussian component is present in a given field. How- 
ever, it does not provide any quantitative information on 
the size of the non-Gaussianity. In this paper, we address 
precisely this issue, by calculating how much the relative 
densities of umbilical points deviates from the universal 
values in the relation to the type and size of the per- 
turbation. In our previous paper [l4j ]. we employed the 
imbalance between the maxima and minima of a field to 
attack the same problem. Besides being applicable even 
when the field itself cannot be observed directly, the ap- 
proach based on umbilics provides an additional probe, 
should the extrema test not be sensitive enough. As an il- 
lustration, consider the case h(r) — H(r)+eH(r) 3 , where 
H(r) is a Gaussian field. Since the perturbation is an odd 
function of H, the symmetry between positive and nega- 
tive values of H is preserved and the relative densities of 
maxima and minima will not differ. By contrast, a study 
of the umbilical points does reveal the non-Gaussianity 
of h, as we will show. 

The outline of this paper is as follows. In sectionlH we 
briefly review the basic properties of Gaussian fields while 
in section UH we introduce the necessary geometric con- 
cepts concerning umbilical points. In sections [Hi] through 
[VT]the various steps and concepts that are needed for this 
calculation are explained, before the final result is arrived 
at in section lVHl The theoretical result is then compared 
to results from computer simulations in section [Villi Fi- 
nally, section ITXl provides a summary and conclusions. 



I. GAUSSIAN FIELDS 

First, we will summarize the main definitions and char- 
acteristics of Gaussian fields that we need; see [l4| for a 
more detailed description. 

A homogeneous and isotropic Gaussian field can be 



defined as 



H{r) = A(k) cos(k ■ 



(1) 



The phases <f>^ are uniformly distributed random vari- 
ables and independent of each other. The amplitude spec- 
trum A(k) depends only on the length of the wave vector 
k and gives the Gaussian field its special characteristics, 
often expressed in terms of the moments of the spectrum 



(2) 



If the spectrum is sufficiently smooth - which we will 
consider to always be the case here - the sum can be 
replaced by an integral, 



K„ 



dfcn(fc)fc™. 



(3) 



Here the integration over the polar coordinate of k has 
been absorbed into n(/c), which is the power spectrum. 
For convenience, we will consider H to be normalized, 
such that K = (H 2 ) = 1. 



II. UMBILICAL POINTS 

Umbilical points are points on a surface where the cur- 
vature of the surface is the same along all directions. The 
curvature depicts how much the surface bends along a 
given direction, just like the second derivative of a one- 
dimensional function does. At an umbilical point then, 
the surface is locally spherical (or flat). 

In order to make a proper mathematical formulation, 
consider a two-dimensional function f(x, y). We consider 
any specific point (xcbZ/o) and any direction given by 
an angle ip. Along this direction, the function can be 
parametrized as 



U(r) = f(x a + r cos -0,2/o + rsin^). 



(4) 



This function now describes what / looks like at (xo,yo) 
along the direction ip. The curvature is the value of the 
second derivative of jf^(r) at r = 0, 



/4'(o) 



d 2 U 



dr 2 



r=0 



= fxx(xc, Vc) cos 2 ip + f yy {x c , y c ) sin 2 ip 

+ 2f xy (x c , y c ) sin tp cos tp 
= (§ + ^cos2ip)f xx + (| - | cos 2ip) fy y + sm2ipf xy 

= \{fxx + fyy) + \{fxx - fyy) COS 2lp + f xy Sm2%(;.(5) 

We can write this in a more lucid form by applying the 
transformation 



\(f xx -fyy)=RcOSa R=^(fxx-fyy) 2 +4f 2 

^fxy 



fxy — Rsina tana = 



fxx fy 



y 

(0) 
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(a) 



(b) 



(c) 



Figure 1: One set of curvature lines (black) around (a) a 
lemon, (b) a monstar and (c) a star. The other set shows the 
same pattern in all cases. The red circle and line segments 
show how the principal direction rotates around the umbilical 
point. 



With this we find 

/#(0) =£(/** + /«,) 

+ ^(f XX -fyy) 2 +^yCOs(2^ - a). 



(7) 



With the curvature now properly defined, we intro- 
duce the two principal directions, which are the direc- 
tions along which the curvature is maximal or minimal. 
The corresponding curvatures are known as the principal 
curvatures. We can easily see from eq. (J7J that these two 
directions are given by 2-0 — a = kn and hence perpen- 
dicular to each other. 

As noted before, at an umbilical point the curvature 
is the same along all directions. In other words, the two 
principal curvatures are the same, and the principal di- 
rections cannot be defined. From eq. ([7)1. the definition 
of an umbilical point is easily seen to be 



fxx — fy 



and 



fxy — 0. 



(8) 



Umbilical points can be classified in three types. The 
distinction can be clearly made when one looks at the 
curvature lines. These are curves which are always tan- 
gent to a principal direction, either the one corresponding 
with the maximal curvature or the minimum one. These 
two sets of curvature lines intersect at right angles, since 
as noted before, the principal directions are always per- 
pendicular to each other. 

At an umbilical point, no principal direction can be 
defined, giving one of the three patterns shown in figured] 
There are three types: lemons, monstars and stars. 

We see that, in each case, the umbilical point is a topo- 
logical defect, having a topological index (see (la ], for ex- 
ample). Formally, the topological index is defined as 



n = — <b Vtj) ■ dl, 
2tt 



(9) 



where the path-integral is taken over an (infinitesimal) 
counterclockwise loop around the defect. In words, it 
counts the number of revolutions the principal direction 
makes when traversing this closed loop. We see from fig.[T] 
that, for the point labeled star, the direction makes half 
a clockwise rotation, which means a topological index 



of —1/2. The minus sign reflects that it rotates in the 
opposite direction with respect to the direction the loop 
is traversed in. The other two umbilical points have index 
+1/2. 

Another characteristic separating the three is the num- 
ber of curvature lines that terminate at the umbilical 
point. For a lemon, this is one, whereas for the other 
two it is three. We see that the third type of umbilical 
point shares properties with both others: it has topologi- 
cal index +1/2, as does a lemon, and three curvature lines 
terminating at it, like a star. This in-between nature of 
the point is reflected in its name: monstar. 

The three types of umbilical point can also be distin- 
guished using the third derivatives, much like the various 
types of critical points can be identified by the second 
derivatives. From eqs. © and ([7]) we see that the two 
principal directions, for which the curvature is maximal 
or minimal, are given by 



tan 2ip 



2 fxy 



fxx f'i 



(10) 



mi 



Note that the directions are given by angles modulo tt, 
so this equation has two solutions: the two principal di- 
rections. The angle 2ip can be pictured as the argument 
of the vector v — ( %f ^ yy ). At an umbilical point, both 
vector components are zero (hence the angle / princi- 
pal direction is not defined). In order to determine the 
topological index, we need to know what the principal 
directions are in close proximity to this point, in order 
to evaluate the infinitesimal loop in eq. ((9]). We can ex- 
pand v using the third derivatives. For a point r near an 
umbilical point ro we have 



fxxx fyyx fxxy fyyy] ( ^ 

^fxyx 2f X yy J \y y u 



A(f- ro). 
(11) 

If A were the identity matrix, then a counterclockwise 
loop around ro would obviously result in 2ip increasing by 
2n, giving index +1/2. In general, A may shear and ro- 
tate f, or may reflect it. The former would have no effect 
on the charge. However, if A includes a reflection, the 
gradient would rotate in the opposite direction and the 
index becomes —1/2. Whether A describes a reflection 
or not is encoded in the sign of its determinant, 

2 det A — (j ^ xxx f xyy) f xyy + {fyyy fxxy) fxxy • (1^) 

Hence, the index of the umbilical point is +1/2 (—1/2) if 
det A is positive (negative). Introducing a = fx X x,P = 
fxxy, J = fx V y,S = fyyy, we thus find (see also [l|) 



a l ~ I 2 + 0$ 




for L, M 
for S 



(13) 



As mentioned before, the criterion separating the 
lemons from the monstars (and stars), is the number of 
(locally straight) lines ending at the umbilical point: one 
for lemons, three for (mon)stars (that is one/three for 
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each principal direction). This can also be expressed in 
terms of a, j3, 7 and 8. Consider again a point r near 
fo- The principal directions are given by "0 modulo ^n, 
hence one of the two is directed toward the umbilical 
point when the argument a of r — fo is equal to ip at r, 
modulo 

To find an algebraic statement of this condition, we 
double both sides: 2a = 2i/j (mod 71"). The right-hand 
side is the argument of v by eq. ((TU)l . We can also find 
a vector whose argument is given by the left-hand side: 
let (y) = f — fo', then 2a is the argument of the vector 

( x 2 ~y ) - this is easily seen by mapping (y) to the com- 
plex number x + iy and taking its square, which doubles 
the argument. The condition for r being on a terminating 
curvature line is that the arguments of these two vectors 
must match modulo n, which translates to I „ y ) and 

' V 2xy > 

A(y) being parallel to each other. This condition can 
be mathematically expressed using the matrix A from 
before and the cross-product, giving 



= | (-2xy x 2 - y 2 ) A 



h-^^)( ( ^;? to ' 3)9 

/3x 3 - {a - 2i)x 2 y + (5 - 2/3)xy 2 - jy 3 . 



(14) 



Note that this equation describes lines passing through 
fo , whereas the curvature lines actually terminate on the 
defect. On one side of fo, this line corresponds with the 
line of maximal curvature, on the other side with min- 
imal curvature. This is easily seen from eq. ([5]), if one 
notes that the second derivatives change sign when pass- 
ing through fo- 

The number of straight lines passing through fo is thus 
equal to the number of (real) roots of this cubic equation 
(that is, by interpreting this as an equation in x/y). This 
is captured by the discriminant: if it is positive, then 
there are three roots; if it is negative, there is only one. 
This results in (see also [l|) 



4(3 7 (a-2 7 )-(<5-2/3) 2 ) 
x (3/3(5 -2/3)- (a - 2 7 ) 2 ) 
-((<5-2/3)(a- 2 7 )- 9/3 7 ) : 

In summary: 



> for M, S 
< for L 
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Umbilical point 


Index +i 


Index — i 


Lemon 


Monstar 


Star 


eq. (fT3|) > 


eq. (jT3| < 


eq. (HI]) < 


eq. CD3) > 



Euler characteristic of the underlying manifold. For the 
two-dimensional plane that we consider, this is simply 
zero. As a consequence, the density of stars (with in- 
dex +1/2) equals the combined density of lemons and 
monstars (with index —1/2). In other words, the star 
fraction, that is the density of stars divided by the to- 
tal density of umbilical points, is always 1/2. There 
is however no topological constraint on the lemon and 
monstar fractions. For isotropic Gaussian random fields 
however, it has been shown that the monstar fraction is 
a M = 1/2 - 1/V5 = 0.053 0, independent of the 
spectrum of the Gaussian field. As a consequence, the 
monstar fraction promises to be a good criterion to test 
the (non-)Gaussianity of a field. Should one be given a 
random field, and find that the monstar fraction is not 
equal to 0.053, one can immediately conclude that the 
field is not an isotropic Gaussian one. 

In the next sections we take a non-Gaussian field h 
that can be described as a Gaussian field H with a per- 
turbation f(H) added to it, and calculate how much the 
monstar fraction aM deviates from the universal value 
0.053 as a function of the perturbation f(H). Our result 
also allows to attack the reverse problem: when given a 
non-Gaussian field of which the type of perturbation is 
known, we can determine the monstar fraction and thus 
reveal the size of the perturbation. 



According to the Poincare-Hopf theorem, for any sur- 
face the total sum of all topological indices equals the 



III. THE GENERATING FUNCTION 



We consider a field of the form h(f) = H(f) + f(H(f)), 
where H(r) is a Gaussian field and / a small nonlinear 
function of H(f) only. As we have seen in eqs. (JSJ, (p~3|) 
and (fT5|) . the monstars can be defined using the second 
and third derivatives of the field h with respect to x and 
y. Determining the monstar fraction thus boils down 
to determining how likely it is that at a specific point 
r the third derivatives a — h xxx (r), /3 = h xxy (f), 7 = 
h xyy (f) and 5 — h yyy (f) are such that eqs. (fT3"| and (fTS")) 
prescribe a monstar, given that the second derivatives 
obey h xx (f) = h yy (r) and h xy (f) = 0. 

In order to determine this, we require the joint prob- 
ability distribution of these seven stochastic variables. 
When we have this, we can set h xx = h yy and h xy = 
and integrate a, /3, 7 and 6 over the appropriate ranges 
to get the density of monstars and all umbilical points 
respectively. The ratio of these then gives the monstar 
fraction. 

We shall arrive at the desired probability distribution 
by determining the corresponding generating function, 
which is defined as the Fourier transform of the proba- 
bility distribution n&i . For a set of n correlated variables 
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{hi} this is 
x(Ai,...,A„) 

= [dh... dh nP (hi, . . . , hn y^ + - +h ^ 



+ q? h n h h h h A ; A /"V. : • ••• (16) 



01,02,33 



Here, the coefficients (. . .) are the moments, or multivari- 
able correlations defined by 

(h n ...h jk ) = J dhx . . .Ah n p(h x , . . . 7 h n )h jl ...hj k . 

Eq. (|16|) is proved by expanding the exponential term by 
term. 

Upon taking the logarithm of x and expanding, the 
quantities known as the cumulants are revealed: 

i 2 

log X = i Ci{hj)\j + Tjy Y C ^ h oi > h h) X ji X 02 



01,02 



(18) 



The cumulants can be written in terms of the moments, 
as can be seen by taking the logarithm of eq. (TIoT) and 
expanding it. For example, 

C{h u h 2 M = (hih 2 h 3 ) - (M(Mi 3 ) - (h2)(h 3 hx) 
-(h 3 }(h 1 h 2 )+2(h 1 )(h 2 )(h 3 ). 

(19) 

In reverse, the moments can be written in terms of the 
cumulants, e.g. 

{hMa) = C(h u h 2 , h 3 ) + C{h x )C{h2,h z ) 

+ C{h 2 )C{h 3 , hi) + C{h 3 )C{h!M) (20) 
+ C(h 1 )C(h 2 )C(h 3 ). 

If all the moments or all the cumulants are known, we 
can construct the generating function and perform an 
inverse Fourier transformation to obtain the probability 
distribution. 

The defining characteristic of Gaussian random vari- 
ables Hi is that all cumulants are zero, with the excep- 
tion of the second-order ones C 2 (Hi, Hj) = (HiHj). In 
this case, the generating function is thus 



(Ai, . . . , A„) = exp CaCHi, -HiJAiAj) . (21) 



The inverse Fourier transformation yields the standard 
distribution for correlated Gaussian random variables 



(see e.g. [13), 
p(H 1 ,...,H n ) 



(27r)™/ 2 v / deto : 



exp(-|^( f r- 1 ) ij ^ j 



(22) 



i,0 



where a is the matrix of correlations, <zy = (HiHj). 

For a Gaussian field H, the derivatives are themselves 
Gaussian fields - therefore the above formula gives their 
joint distribution. For the non-Gaussian field h, there 
are some small corrections to this distribution. To find 
these corrections to first order, we need to determine the 
cumulants to first order in f(H). We will see that only 
a small number of cumulants are nonzero up to this or- 
der. Before we proceed to derive them, we switch to a 
complex coordinate system which allows for optimal us- 
age of translational and rotational symmetry, which h 
has inherited from H for the type of perturbations under 
consideration. 



IV. COMPLEX COORDINATES 
REPRESENTATION 

To find the distribution of umbilical points, we now 
have to find the joint distribution of the seven second 
and third derivatives of h. All these variables can be 
combined into a more compact form by using complex 
coordinates. These will make it easier to evaluate the 
integral that determines the monstar density, and will 
help us to work out the probability distribution with the 
help of symmetry. 

The complex coordinates are given by 



z = x + iy 
z* = x — iy 



x = \(z + z*) 

y = \i{z* - z) (23) 



Of course, as complex numbers, z and z* are not indepen- 
dent; however we can formally define partial derivatives 
with respect to each of them, using the chain rule, just 
like we could if this transformation involved a real num- 
ber instead of i. 

The derivatives with respect to z and z* are given by 

A = iA_i l A i A = iA + iiA. (24) 

dz 2 dx 2 dy ' dz* 2 dx 2 dy 

We see that the derivatives with respect to z and z* are 
each other's conjugate, but again, we consider both to be 
linear transformations of d x and d y . 

The usefulness of using z and z* can be immediately 
seen from h zz : 

h zz = d z h = j(d x -id y ) 2 h = j(h xx - h yy + 2ih xy ). (25) 

We see that the definition of an umbilical point can be 
captured in one equation: h zz — 0. 
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The various types of umbilical points were defined in 
eqs. (fT5|) and (fT5)) using the "normal" third derivatives 



hzZZi hzZZ* 7 h<ZZ*Z* 

monstar become 



P; h x yy 

and h. 



7 and h y 



5. In terms of 



the two conditions for a 



\h Z zz*\ 2 -\h Z zz\ 2 >0, (26a) 
27\h zzz \ 4 - \h zzz , | 4 - l%\h zzz \ 2 \h zzz , | 2 

- A(h zzz h z zz , z , + h x . x . x .h*„.) > 0. (26b) 



represent h zzz h z * z * z * and 



Here \h zzz \ 2 and \h zzz , 
h zzz *h zz , z * respectively. 

The density of monstars will essentially be given 
by integrating the probability distribution p(h zz = 
0, h ZZZl h zzz * ) over the range defined by these conditions. 
This probability distribution is determined by the cumu- 
lants of combinations of the three variables. Rotational 
and translational symmetry however imply that only a 
small number of these combinations yield a nonzero cu- 
mulant. 

First, consider the consequences of the isotropy (ro- 
tational symmetry) of the field h(f) for a moment like 
(h z *(r)h zz (f)). Note that, due to homogeneity (trans- 
lational symmetry), this moment does not depend on r; 
it will often be dropped from now on. Isotropy implies 
that this moment should not change if we rotate the field 
around r, over any angle a. In terms of z and z* , this 
results in the transformation 



d z , = 
cU = e la d z 



(27) 



As a result, we get (h z >h z '* z >*) = e la (h z h z * z *). Since we 
argued that the two expectation values must be equal, 
for any a, we must have {h z »h zz ) = 0. 

In general, following a rotation expectation values pick 
up a factor e lka , where k is the number of z* deriva- 
tives inside the bracket minus the number of z deriva- 
tives. By the above argument, the expectation value is 
zero if fc 7^ 0. Therefore, an expectation value can only 
be nonzero if the numbers of z and z* derivatives inside 
the bracket are equal. Since a cumulant is a sum of prod- 
ucts of expectation values, featuring every variable once 
in every product (compare eq. (I19|) V the same property 
applies to cumulants. 

The homogeneity (translational symmetry) of the fields 
under consideration provides another useful trick that 
relates different cumulants to one another. As already 
stated, a moment like (hi(r) . . . h n (r)) does not depend 
on r. Hence the derivative of this with respect to z or z* 
is zero. Applying the product rule: 



Q = d z {hx...h n ) 
= ((d z h 1 )h 2 ■ ■ ■ K) + . . . + {hih 2 ■ ■ ■ {d z h n )). 

For n — 2, this gives the useful relation 

((cUi)M = -{hi{d z h 2 )). 



(28) 



In essence, for a two-point correlation it is possible to 
"transfer" a z derivative from the one term to the other 
at the cost of an overall minus sign. The same applies 
of course to a z* derivative. For example, we find the 
relation (h zz h z * z *) = ~{h z h zz * z .). 

Together, these two symmetries constrain the proba- 
bility distribution p(h zz , h ZZZl h zzz »). In particular, they 
explain why the monstar fraction is always the same for 
any Gaussian distribution [lj. A Gaussian distribution 
does not have many degrees of freedom to start with; 
only the two-point correlations between the variables are 
adjustable. In this case, the two-point correlations be- 
tween any two of these variables is zero, by rotational 
symmetry, while the variances of h zzz and h zzz * are 
equal by translational symmetry. Hence (after setting 
h zz = to identify the umbilical points) , the distribution 
p(h zz = 0,h ZZZl h zzz *) is always the same apart from a 
scale, and that determines the monstar fraction. This 
argument can be generalized to singularities in the po- 
larization field of light (even though the field might not 
be derived from a scalar field h), and so Gaussian polar- 
ization fields have the same monstar fraction as well, as 
shown in Q. 

On the other hand, when h has non-Gaussian contri- 
butions, there are many more cumulants, and symmetry 
is not enough to constrain them any more. For the field 
h = H + f{H) we are studying, we proceed to calculate 
the cumulants explicitly. 



V. THE CUMULANTS 

Although the problem is now cast in terms of complex 
derivatives, the recipe outlined in section Hill still applies. 
The task is to determine the cumulants of h zz , h ZZZl h zzz * 
and their conjugates up to first order in the perturbation 

/• 

These cumulants have the general form 
C n {D\h, . . . , D n h), where each Dj represents a num- 
ber of z and z* derivatives. For the moment, let us 
consider each Djh to be at a different point fj, i.e. 
C n (Dih(fl), . . . , D n h(fn)). Later we will set all points 
equal again. For convenience, we shall drop the vector 
notation, i.e. r» = fl. Since now each derivative Dj acts 
only at a specific point, we can bring them outside the 
cumulant: 



C n (D 1 h,...,D n h) 

= D 1 ...D n C n (h(r 1 ), 



, h{r n )) 



(30) 



Let us write h(rj) = hj for shortness, and focus on 



C n {h x 



Inserting hj = Hj + f(Hj), expanding 



(29) the cumulant and keeping only terms up to first order in 



7 



/ yields 

C n {h\, . . . , h n ) = C n (Hi, . . . , H n ) 

+ Cn(f(Hi),H2, ■ • ■ , H n ) 
+ C n (H 1 J(H 2 ),...,H n ) 
+ C n {H\,H2, ■ • ■ , /(.Hn)). 



(31) 



The first term on the right-hand side is now simply the 
cumulant of a set of Gaussian random variables, which, as 
discussed before, is zero for n > 2. The other terms can 
be evaluated perturbatively and are equivalent to each 
other. Consider the second term as an example. For a 
cumulant involving Gaussian variables and one function 
of a Gaussian we have (see Appendix EJ: 



Cn(f\Hi),H2, • ■ • , H n 



{f^- 1 \H 1 )){H 1 H 2 ){H 1 H z ) . . . (M„>. 



(32) 



When we reinsert the derivatives D2 through D n from 
eq. (|30|) and set r 2 = . ■ ■ = r n — r we get 



n—r 
(33) 



D 1 ...D n C n (f(H(r 1 )),...,H(r n )) 

n — ...—r n 

= D 1 {f^- x \H 1 )){H l D 2 H) . . . (HiD n H) 



Now we can reinsert D\ and then set r% — r, as pre- 
scribed by eq. (|30|) . Remember that D\ only acts on H\ . 
Due to the product rule, we have to consider all possible 
ways in which the derivatives in D\ can be distributed 
over all Hi's. Recall from section HVl that, after setting 
7*i = r, each expectation value can only be nonzero if the 
number of z and z* derivatives inside are equal. Note 
also that (f^ 1 ^ (Hi)) does not depend on n, hence any 
derivative of it is zero. Therefore, the only nonzero con- 
tributions stemming from the product rule are those dis- 
tributions that make the number of z and z* derivatives 
equal inside each bracket. 

We consider the cumulant C^(h ZZl h zz * z *, h zz * z ») as an 
example to demonstrate the procedure. First, we find 

Cz(h zz , h zz * z * , h zz * z *) 

— C%(H ZZ , H zz « z », H zz » z *) 

+ d ZlZl ((f"(H))(H 1 H zz * z .)(H 1 H zz * z .)) (34) 
+ 2d Z2Z , z; ((f"(H))(H 2 H zz )(H 2 H zz * z *)). 

The first term is zero, since all cumulants of Gaussian 
variables are zero beyond second order. For the second 
term, we need to consider how to distribute the two d Zl 
derivatives to make all expectation values nonzero. The 
only possibility is to put one d Zl in front of each H\. 
Note however that this term appears twice in the product 
rule, because there are two ways of distributing the two 
derivatives. After setting r\ = r we thus have 



d ZiZi ((f"(H)) (HxH zg * z .)(HxH,*.*.)) 
= 2(f"(H))(H z H zz * z .) 2 . 



In the third term on the right-hand side of eq. ([33)1 we 
have one d z and two c^.'s to distribute. The first H 2 
needs <9 Z * Z * to balance the derivatives and the other takes 
the d z derivative. There are no multiple ways to dis- 
tribute these derivatives in this case and therefore 

d Z2Ziz ' 2 ((f"(H))(H 2 H zz )(H 2 H zz , z *)) 

= (f"(H))(H z , z ,H zz )(H z H zz , z ,). (36) 

Combining everything together results in 



Cs(h zz , h zz * z * , h zz * z *) 

= 2(f"{H))(H z H zz , z ,){(H z H z 



\H Z * Z *H ZZ )). 

(37) 



Finally, due to translational symmetry, we have 
(H Z H ZZ , Z ,) + (H Z , Z ,H ZZ ) = d z (H z H z , z ,) = 0. We thus 
find that C 3 (h zx ,h zz * z *,h xz * z *) = 0. 

Now we will show that there are only a finite number of 
nonzero cumulants (up to first order in /). In fact, there 
are none beyond the fourth order. Consider eq. (|31[) with 
n > 4. The first term (zero order) is zero because it is 
the cumulant of more than two Gaussian variables. For 
the other ones we apply the recipe of eq. (|3"3"1) . We have 
n — 1 brackets in which the z and z* derivatives need to 
be matched. Since we are only considering the variables 
h zz , h zzz , h zzz * and their conjugates, each one has a 
mismatch to begin with. However, since D\ has only 
three derivatives at most, it is not possible to balance 
the derivatives in all n — 1 brackets. 

This "lack of derivatives" also kills a lot of cumulants 
of lower order, especially fourth order. For example, 
in C4(h zz ,h z * z *,h zzz *,h zz * z *) the first two variables re- 
quire two derivatives to balance the derivatives and the 
other two require one. Therefore, no matter from which 
variable the derivatives are distributed, there is always a 
shortage. 

All the nonzero cumulants are listed in table |Vl Two 
parameters were introduced: 



o" = (H ZZ H Z , Z *) — —(H Z H ZZ * Z *), 
t = (H ZZZ H z * z * z *) — (H zzz * H zz * z * 



(38a) 
(38b) 



(35) 



Note that the trick based on translational symmetry was 
used to equate the expectation values. In the second 
equation, it was used twice (transferring a d z one way 
and a d z * the other way). 

The parameters a and r are related to the moments 
K n of H . This is most easily accomplished by writing H 
in complex variables: 

H = ^2 A(k) cos(fc • r + ^) 
% 

= J2 MM) cos (\{k*z + kz*) + (j> k ). (39) 
fc 

Here k = k x + ik y is the complex analogue of fc = (^ x ). 
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Ci{h 7 . 7 ., h z * z *) 

C^ihzzzi h z *z*z*) 
C2(hzzz* , hzz'z' ) 

C 3 (hzz, hzzz*, hz'z'z*) + conj. 

C^ih Z zz* f h Z zz* i h Z z*z* •> h Z z*z* ) 

C^{h zzz i h zz * z* 5 h zz * z * , h Z z* z* ) ~h conj. 



<r(l + 2{f'(H))) 
r(l + 2(/'(fl)» 
r(l + 2(/'(fl)» 
-3a 2 (/"(//)) 
-8a 3 {f'"(H)) 
-6a 3 {f"'(H)) 



Table I: All nonzero cumulants. The two asymmetric cumu- 
lants have a conjugate twin in which all z's and z*'s are in- 
terchanged. 



With this we find 

H zzz = J2 A (k)l(k*) 3 sin {\{k* z + kz*) + fa), 

k 

(40a) 

H z . z * z . = J2 A (k)lk 3 sin {\{k*z + kz*) + fa). (40b) 

k 

Hence: 

r = (H ZZZ H Z , Z , Z ,) = (i(fc*) 3 |fc 3 ) = ±K 6 . (41) 
Similarly, we have a = j§K^. 



VI. PROBABILITY DISTRIBUTION 

With the aid of the cumulants we can build the loga- 
rithm of the generating function (see eq. (|18|) ) . provided 
that we identify the appropriate variables in Fourier 
space. Consider h zz and h z * z * for example. These com- 
plex variables represent two real variables £ x and £ y , the 
real and imaginary part of h zz . Let A^ and X y be their 
Fourier counterparts. The generating function is 

X(K, A,,....; / d^ y . . . £„• • . .)e^ A *+«^+-> 



The exponent can be written in terms of the complex 
variables, 



£xX x + £y^y — hzz^zz + h z * z *X z 



(43) 



where we define X zz = ^(A^ + iX y ). Then X zz is the 
complex Fourier variable corresponding to h z * z * and we 
likewise introduce X zzz and X zzz * , which are conjugate 
to h z * z * z * and h zz * z *. We will define integrals with re- 
spect to the complex Fourier variables, e.g. with respect 
to d 2 h zz , as integrals over the real and imaginary parts of 
h zz , and the inverse Fourier transform will be performed 
by integrating over the real and imaginary parts of the 
A's. 



The generating function is thus 
logx = - C 2 (h zz ,h z * z *)X z * z *X zz 

C*2 {hzzzi h z * z* z*')Az* z* z* ^zzz 
C*2 {hzzz* i h Z z* z* ~)^zz* z* ^zzz* 

iC*3 {fizz-; hzZZ* 7 h z * z * Z* ) ^Z* Z* ^ZZ* Z* ^ZZZ 

— iC^{hzzi h Z zz* j h z * z * Z *)A ZZ \ ZZZ * \ z * z * z * 

~\~ "^C^hzzz* 1 h Z zz* 1 h Z z* z* 1 hzz* z*~)^zz* z* ^zzz* 
g {h Z zz ; h Z z* z* 1 h Z z* z* 1 h Z z* z*~)^z* z* z* ^zzz* 

~t~ Q C4 {hzZZ ; h Z Z* Z* 1 h Z Z* Z* 1 h Z Z* Z*^)^ZZZ ^ZZ* Z* ' 

(44) 

Upon entering the cumulants from table [Vj 

1°SX = — &X ZZ X Z » z * — t(X zzz X z * z * z * + X zzz * X zz * z *) 
+ ZicT 2 {f"(H))(X zz X zzz »X z , z * z » 

+A Z » 2 * X zz * z » X zzz ) (45) 
- ^ (f"'(H))(2X zzz , X 2 ZZ , Z , + X ZZZ X\ Z , Z , 
+A Z * Z » Z *A ZZZ ,). 

Here a = a(l + 2{f{H))) and f = r(l + 2(f'(H))) have 
been introduced. The factors in front of the cumulants 
are the factor i k jk\ in eq. ()18[) multiplied with the number 
of permutations of the A's. 

To obtain the probability distribution, we take the ex- 
ponential and perform the inverse Fourier transformation 
(see eq. (p~6|) ) . This gives an integral over the exponen- 
tial of a polynomial of degree 4. However, all terms of 
degree three and four are of order /, so we can expand 
the exponential and be left with only square terms in the 
exponent. The result is 

X = exp ( — bX zz X z * z * — f(X zzz X z * z * z * + X zzz * X zz * z * )) 
x (l + 3ia 2 (f"(H))(X zz X zzz *X z « z * z * 

U X z * z * X zz * z * X zzz ) 

- 0-3 (/"' \H))(2X 2 ZZZ »X 2 ZZ , z , + X ZZZ X\ Z » Z * 

+W,.AL,.))- (46) 

Now we can take the inverse Fourier transform. Note 
that X zz and A z * z * are each other's conjugate. Upon inte- 
grating the real and imaginary parts of X zz and imposing 
A z * z * = A zz (the same procedure applies to the other two 
pairs of A's), one obtains 

p{h zz t h z * z * 1 h zzz , h z * z * z * , h zzz * , h zz * z * ) 

f d 2 X zz d 2 X zzz d 2 X zzz * 
Zr X 



xe 



—i(X zz h zz +X zzz h zzz -{-\ zzz *h zzz *-\-cony) 



• (47) 



Note that the denominator is 7r 6 rather than (2tt) 6 be- 
cause of the factor of h in the definitions of the A's (see 
eq. (O). 

The Fourier transform of a Gaussian function multi- 
plied with a polynomial is easy to perform by noting that 
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multiplying by A in Fourier space is equivalent to taking 
a derivative in normal space: 



dAA"/(A)e 



-iXh 



t|-) /<1A/(A),-' A ''. UN) 



The inverse Fourier transform of the Gaussian part of the 
generating function is 



12 1 ^ rd 2 A 



f d 2 X zz d 2 X z 

J IT' 



,~~ 2 exp (-\\h zz \ 2 - 



? + \h A 2 ' 

ZZZ | ~ \ IL ZZZ | , 



(49) 



and the final result reads 

pifozz i h z * z * 3 h zzz j h z * z * z * , h zzz * , h zz * z * ) 

1 - 3a 2 (f"(H))^ (2$t(h z . z .h zz . z .h zzz )) 



2|/ w | 4 + 23?(/i zzz /4, z .) 



1 / 
q--2 exp I 



— )-(50) 



VII. MONSTAR FRACTION 

Once the joint probability distribution of the relevant 
derivatives is obtained, we can set h zz = h z * z * = 0, which 
defines an umbilical point. The joint probability distri- 
bution states how likely it is that h zz and h z * z * are close 
to zero for a certain point r. What we need however is 
for h zz and h z * z * to be exactly zero for a point close to 
r, since we are looking for a density with respect to the 
(x, y)-plane. For this, we need to go from a probability 
density with respect to h zz and h z * z * to one with respect 
to z and z*. This is accomplished by multiplying p with 
the Jacobian 



J 



d{h zz ,h z * z *) 



d(z,z*) 



\h z 



(51) 



The last step is to integrate this product over h zzz and 
h zzz * , either over all possible values, or just over those 
satisfying eq. (|2"rj|) , to get the density of umbilical points 
and the density of monstars respectively: 



/ d 2 h zzz d 2 h zzz , p(h zz =Q,h zzz ,h zzz *) 

JR. 



(52) 



x J(h zzz , h zzz *) , 



where R represents the range of integration: the entire 
space to get the density of all umbilical points, or eq. (|26|) 
for just the monstars. 

First we simplify by introducing polar coordinates, 



J<t> 



h z 



\K 



(53) 



Next, we introduce 

\h zz 

U = T- 
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(54) 



We find that we can rewrite the two conditions for mon- 
stars, eq. p6p. in terms of u and 5 only: the first one is 
simply u > 1, while the other is 

< 27 - u 4 - 18u 2 - 8u 3 cos 25 
= (3 - u) 3 (l +u)- 16u 3 cos 2 S 



cos S < 



16u 3 



(55) 



Since the fraction on the right-hand side is negative for 
u > 3, we can extend the first condition to 1 < u < 3. 

The fact that the monstar conditions depend only on 
u and S can be understood as follows: the type of um- 
bilic should not be affected by rescaling and/or rotating 
the plane. Rescaling would add the same (real) factor 
to h zzz and h zzz », hence the type of umbilic should, as 
far as the moduli are concerned, depend only on the ra- 
tio \h zzz \/\h zzz *\. A rotation introduces phase factors as 
given by eq. (f27|) . We see that a rotation over an angle a 
causes h zzz to pick up a factor e~ 3ta while h zzz * picks up 
e~ la . Therefore, the only combination of <j> and 9 that is 
invariant under rotations is 39 — <fi. 

Now we return our attention to the probability distri- 
bution. First we rescale h zzz and h zzz * , 



(56) 



This leads to 



p(h zz = 0,v,w)J 



:(l - ^{f"'(H))(A - 8M 2 + 2M 4 + vw* 3 + v*w 3 )^j 



xe 



Here we dropped an overall coefficient, which is of no 
importance since we are only interested in the ratio of 
the densities of monstars and all umbilical points. Note 
that f now only appears in the term proportional to /. 
Since we are not interested in higher orders of /, we need 
only consider the leading order of f, which is r. For 
convenience, let us define 



e=" </'"(ff)>. 



(57) 



Furthermore, note that multiplying p with the constant 
1 + 4e - which we may do since we are only interested in 
the density of the ratios - causes the 4 inside the paren- 
theses to be canceled out (up to first order). 

Next, we move to polar coordinates, as we did before: 

m 



pe 



itj) 



W 



re 



(58) 
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Figure 2: The monstar fraction o.m of H + eH 3 as a function 
of e, where H has a disk spectrum (/x = ||). The data points 
stem from simulations, the solid line is eq. (|75|l . 



and then substitute r — up and (f> 
transformations we have 



3(9 - 25. With these 



n oc / p u dpdud8dS e 

„2„.2 



" ( " +1) p 2 \u 2 



x (l - e(-8pV +2p 4 (u 4 + u 3 cos 25)). (59) 

Finally, we integrate over p and 9 to find the probability 
distribution p(u, 5). The integration over 9 simply gives 
a factor of 2ir, while the integral over p has the form of 
a polynomial times a Gaussian. For this we can use 



dpp 2n+1 e-P (" +1 > = 



o 



The result is 



p(u, S) oc 



II 



1 + 24 



2(u 2 + 1)™ +1 ' 
_u 2 (l — u cos 26) 



(60) 



(61) 



(w 2 + l) 3 V (u 2 + I) 2 

The monstar density is proportional to the integral of 
p(u, S) over the range 



1 1 -l / / \° ~ "rl 1 t u \ r 1 
1< M <3, cos 1 ^ ^ '-)<6<\n, 

(62) 

while the total density of umbilical points is proportional 
(with the same prefactor) to the integral over the range 
< u < oo, < 8 < |7r (extending the integration 
range of S from ^tt to 2ir would just add a factor of 4 to 
both integrals). The latter can be done analytically: the 
integration over d is trivial, while the remaining integral 
over p can be split into two parts which, apart from a 
factor sgn(w 2 — 1), are both of the form 



du 



u 2 -1 



(u 2 + l) 2 \u 2 + 1 



n + 1 \ u 2 + 1 



n+l 



(63) 



Figure 3: The monstar fraction qm of H + eH 3 as a function 
of e, where H has a Gaussian spectrum (u = |). The data 
points stem from simulations, the solid line is eq. (|75|l . 

With this we find 



ntot oc ~tt 



,4 1 1 



24e-J 



2 (u 2 + l) 2 4 (u 2 + l) 4 

24e- 1 



2 (u 2 + l) 2 4 (u 2 + l) 4 



§*r(l + 3e). 



(64) 



For the monstar range, the integral over S can be per- 
formed. The integral over the cosine gives 



dS cos 26 — 



1 



16u 3 



/ cos 1 (. . .) 

All together: 

r\ u(u 2 - 1) . _j 

n M oc / du sin 
24u 2 



^(u 2 - l)(9-u 2 ) 3 . (65) 



(u 2 + l) 2 



; (3-m) 3 (1+m) 
16u 3 



'(3-m) 3 (1 + m) 
16m 3 



3 v /(m 2 -1)(9-u 2 ) 3 
2(u 2 + l) 2 



= h+eh- 



(66) 



The integrals can be done numerically. The monstar 
fraction is then 



o>m 



n M 



h + eh 



ntot g7r(l + 3e) 7r 
0.053 + 0A29fi{f"'(H)) + 0{f 2 ) 



- = ^-{h + e{I 2 - Zh)) + 0{e 2 ) 

(67) 
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where 



^^~2=J^ (0<M<1) 



(68) 



Note that the zeroth order result matches the one in [l[ . 

Remember that we set Kq = (H 2 ) = 1 for convenience; 
if we drop this condition, then Kq enters the denominator 
of the expression above. 

There is an alternative expression for the term 
(f"'(H)) in eq. (jBT)) . Since H is Gaussian, with mean 
and deviation Kq — 1, we can write 



(f"(H)) = \&zf'\z)e-^l\ 



(69) 



Repeated partial integration yields 

(f"(H))= [dzzf"(z)e~* 2 / 2 

dz (z 2 - l)f'{z)e- z2 ' 2 



= ydz (z 3 - 3z)/(z)c 
= {(H s -3H)f(H)). 



-z 2 /2 



(70) 



The kurtosis of a stochastic variable is defined as the 
fourth cumulant divided by the square of the second, 
which gives 



- (h 2 ) 2 6 - 

If we enter h = H + f(H), we find 



(71) 



(if 4 



4<if 3 /(ff)) 



(H 2 )+A(Hf(H)) 
_ 4(if 3 /(if)) - 12(Hf(H)) 

l+A(Hf(H)) 
= 4(H 3 f(H))-l2(Hf(H)). 



3 + 0(f) 
0(f) 



(72) 

We see that (f"'(H)) = k/4 up to first order, which 
remains true if Kq =/= 1. Hence an alternative form of 
eq. (|6T)) is 



a M = 0.053 + 0.107/iK + 0(f), 



(73) 



where k is the kurtosis of h. 

By comparing the fraction of monstars in a given field 
to the formula just found, we can determine one pa- 
rameter of the deviation from a Gaussian distribution, 
(f"'(H)). This assumes that the field h is given by 
h = H + f(H). To test this, one could, if possible, also 
measure the distribution p(u, S) to test that it has the 
right form, eq. (|6T|) . Measuring p(S) only could also suf- 
fice. For a Gaussian field if, all values of 5 should be 
equally likely, whereas integrating eq. (1611) shows that 
the distribution we expect for h is 



1 



p(S) = -(l-4£cos2<5), 
where we define 5 to lie between — ir/2 and ir/2. 



(74) 



VIII. COMPARISON WITH SIMULATIONS 

The most basic example of a non-Gaussian variable for 
which cq. (|67|) can be tested is h = H + eff 3 , for which 
(f"'(H)) = 6e. Then we have 



a M = 0.053 + 2.576^£ + 0(e 2 ). 



(75) 



Eq. ([57)1 was compared to results from simulations. Gaus- 
sian fields if of square size were generated by adding to- 
gether a large number (a few hundred) of waves with ran- 
dom phases, in the spirit of eq. ([1]). We then counted the 
number of monstars and umbilics in h = H + eH 3 for var- 
ious values of e. Periodic boundary conditions were ap- 
plied to reduce finite-size effects. Furthermore, we chose 
spectra for which the spectrum decays very quickly (or is 
zero) for large k: a disk spectrum, 



A(k) 2 ~ 9(k - k ) K 2n = 



Kq 



(76) 



and a Gaussian spectrum, 



A(k) 2 - exp(-/c 2 /2fc2) 



K 2n = 2 n n\kl n 



(77) 

The simulations were done in the same way as described 
in Q. 

A very good agreement between theory and simulation 
was found for both spectra (see figs. [5] and [3]), for e up to 
about 0.01. For larger values of e, nonlinear terms start 
to dominate. 

Another thing to note is the sensitivity: in eq. (|75[) we 
see that the prefactor of the perturbation term is very 
large compared to the leading order. As a result, even 
for small e, the relative deviation from the universal 0.053 
is quite large, as can be seen in the graphs. Therefore, 
measuring the monstar fraction of a given field proves to 
be a good method for detecting and quantifying small 
deviations from Gaussianity. 



IX. CONCLUSION 

We have calculated how the density of monstars 
changes for a non-Gaussian field given by h = H + f(H), 
where if is a Gaussian field. Comparing our formula to 
data allows to measure the parameter (f"'(H)) = k/4 
of the non-Gaussian contribution. Furthermore, one can 
also measure the distribution of the parameters u and 
5 that define the type of umbilical point. The expected 
distribution is eq. (I5IT) . or eq. (|74l) for just the variable 
8. Measuring these distributions further constrains the 
value of (f"'(H)), and more importantly, it gives a test 
that the non-Gaussianity really arises from a local non- 
linear transformation of a Gaussian field. 

Even though in general the derivatives of h have an 
infinite number of nonzero cumulants up to first order 
in /, it turns out that the cumulants of the variables 
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which are relevant for the umbilics (h zz , h zzz and h zzz * 
at a single point) vanish beyond the fourth order due to 
symmetry. As a result, we found the interesting result 
that (up to first order) olm depends only on (f"'(H)). 

For a more general type of non-Gaussian field there 
would be more independent variables and hence more 
nonzero cumulants. However, it often still holds that the 
higher order cumulants are of less importance. In this 
case, one can consider only the cumulants up to a specific 
order (e.g. fourth order), of which still many would be 
zero due to symmetry. Applying the same procedure as 
outlined here could then reveal the monstar fraction up 
to first order. 
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Appendix A: Proof of eq. ({32]) 

In this section we prove the identity 

C n (f(Hx),H 2 ,...,H n ) 

= (f^- 1 HH 1 ))(H 1 H 2 )(H 1 H 3 )...(H 1 H n ), 

for Gaussian variables Hi. 

Recall the definition of a cumulant, eq. (fT5)). The gen- 
erating function \ is the Fourier transform of the proba- 
bility distribution (see eq. ()16[) ). which for Gaussian vari- 
ables is eq. (|22|) . This leads to the identity 



C n (f(Hi), H 2 , . . . , H n ) 



expC-lSijgj/frjhj) 
(2tt)™/ Vdct cr 



(A2) 



Ai = 



=A„=0 



First, the integration over h 2 through h n is performed. 
This partial Fourier transform is not trivial. If h\ were 
included, the answer would be simply eq. (|2T|) . We can 
however use this result and take the inverse Fourier trans- 
form of it with respect to Ai to get the desired result: 



dh 2 ...dh n e *(^+...+A 



(27r)"/ 2 Vdeto : 



ij A j Xj 



dA_i 
"2T 



exp 



: exp 



i>2 

(ihi +^2a lj X j y 

i>2 



2cti 



3 X! rr '-' A ' A 

i,j>2 



(A3) 



The integration was performed by completing the square. 
Now we do the Fourier transform with respect to hi and 



include the logarithm present in eq. (IA2|) . which leads to 



log J dhl...dh n e ^f{hi)+X2h 2 + ...+X n h n ) 

exp(-^E IJ ^ 1 ^fej) 
(2?r) n / 2 Vdet a 

= log /d/jie iAl/(,ll) — L 



e 2^jl ih i+^-.j>2 CTi j a j'] 2 (A4) 



i,3>2 



In accordance with eq. (|A2[) , we must take the derivative 
of this equation with respect to the A's and set them to 
zero. First the derivative with respect to Ai is taken. 
This causes the first term to vanish, since it does not 
depend on Ai. This simplification can be regarded as the 
main reason why the final result depends on / in a rather 



13 



simple way. What remains is 



This results in 



dXi J V27rcrn 

e 2^j[ i?1 l+E J >2 <J lJ A j] 2 



Ai=0 



1 



d^i/(/ii)e^ [l ' ll+E ^ 2<TljA312 . (A5) 



V27rcrii 

For each with k > 2 the derivative yields 



■• 9 p2^j['fel+Ej>2^^1 2 



-i u 



/i 1 +^<Ti J -A j )e^r [ihl+E ^ 2<Tl3 ^ ]2 



J>2 



(A6) 



Applying this for all k and subsequently setting all A^ to 
zero then gives 



1 j d\ 2 '"dx n 



CTll 



n-1 



n 1 ' 

J> 2 _ 



-h 2 1 /(2a 11 ) 



(A7) 



C n {f{H 1 ),H 2 ,...,H n ) = I 

S>2 



(A8) 



Integrating by parts n — 1 times leads to 



c n (f(Hi),H 2 ,...,H n )= ( J : 



(A9) 



V27TO-11 7 



Finally, we identify the integral (along with the prefactor) 
as the expectation value of and try = (HiHj), 

which gives us 



Cn(f(Hi), H 2 , H n ) 

= {f^- 1 \H 1 ))(H 1 H 2 ){H 1 H 3 ) . . . (HiH n ), 

the equation we set out to prove. 



(A10) 
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